Counterion-mediated Electrostatic Interactions between Helical Molecules 
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We study the interaction of two cylinders with helical charge distribution mediated by neutral- 
izing counterions, by analyzing the separation as well as the azimuthal angle dependence of the 
interaction force in the weak and strong coupling limit. While the azimuthal dependence of the 
interaction in the weak coupling limit is overall small and mostly negligible, the strong coupling 
limit leads to qualitatively new features of the interaction, among others also to an orientationally 
dependent optimal configuration that is driven by angular dependence of the correlation attraction. 
We investigate the properties of this azimuthal ordering in detail and compare it to existing results. 
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I. INTRODUCTION 

Electrostatic interactions often play a dominant role 
in biological and soft-matter systems [l|. In aqueous 
environments charges on macromolecular surfaces such 
as membranes, self-assembled micelles, globular proteins 
and fibrous polysaccharides tend to dissociate and af- 
fect a wealth of functional, structural and dynamical 
properties [2!]. Charged biological macromolecules and 
self- assembled colloidal systems also show a ubiquitous 
patterning of the structural charges that are seldom dis- 
tributed uniformly on the macromolecular surfaces. This 
can be most clearly seen in the patchy nature of charge 
distribution on various proteins, the notorious double- 
helical charge motif seen on DNA [2] and then all the 
way to self-assembled mixed charged bilayers [4] and dec- 
orated carbon nanotubes (CNT). In the latter case too, as 
was realized recently @, 0] j one can have a pronounced 
helical motif of the charges when CNTs associate with 
single stranded DNAs that charges it up and makes this 
quintessential hydrophobic macromolecule soluble in wa- 
ter. 

Interactions between macromolecules exhibiting helical 
charge motifs were investigated vigorously by Kornyshev 
and Leikin ^ who formulated a Debye-Hiickel-Bjerrum 
type theory for interactions between DNAs. In their ap- 
proach the charge distribution on DNAs, including the 
counterions in close proximity to the charged phosphate 
groups, is assumed a priori and the interaction between 
two such molecules in close proximity is calculated within 
the Debye - Hiickel approximation based on this assump- 
tion. In the case of ss-DNA covered CNTs there were 
various attempts to derive the presumed helical charge 
motif based on either approximate Debye - Hiickel (DH) 
arguments or MD simulations 0]. Here the interac- 
tions between decorated CNTs are less well studied since 
it is not clear whether the helical charge motif would 
not be significantly perturbed by the interaction itself. 
However, in both cases of macromolecules with helical 
charge motifs, the electrostatic interaction is basically 
treated on the DH level and should thus be valid asymp- 
totically for sufficiently large separations between macro- 



molecules |10l |. There are other regions in the parameter 
space of these type of systems where one expects 
the DH framework to break down and one should see the 
emergence of completely different type of physics. Unfor- 
tunately we are still not close to having a consistent and 
possibly asymptotically exact formulation of the statisti- 
cal mechanics of charged systems, containing a mixture 
of salts, polyvalent counterions and strong fixed charges, 
though there is some recent progress in this direction [lO[ . 

In the absence of a general approach that would cover 
thoroughly all the regions of the parameter space one 
has to take recourse to various partial formulations that 
take into account only this or that facet of the problem. 
In this respect the counterion-only or the one-component 
Coulomb fluid model system has proved to be of substan- 
tial value [ll|. Heuristically as well as numerically. A 
proper understanding of the behavior of charged systems 
would thus start with the analysis of counter-ion distribu- 
tion around charged macromolecular surfaces, neglecting 
completely the effects of salt. The traditional approach to 
these one-component Coulomb fluids has been the mean- 
field Poisson-Boltzmann (PB) formalism applicable at 
weak surface charges, low counter-ion valency and high 
temperature [l3,[ili[3|- The limitations of this approach 
become practically important in highly-charged systems 
where counterion-mediated interactions between charged 
bodies start to deviate substantially from the mean-field 
accepted wisdom [ll|. One of the most important re- 
cent advances in this field has been the systematization 
of these non-PB effects based on the notions of weak 
and strong coupling approximations. This approach has 
been pioneered by Rouzina and Bloomfield [l^ and elab- 
orated later by Shklovskii et al. [l6| . Levin et al. 
and brought into final form by Netz et al. [ill . \\A 



These two approximations allow for an explicit and ex- 
act treatment of charged systems at two disjoint limiting 
conditions whereas the parameter space in between can 
be analysed only approximately and is mostly accessible 
solely via computer simulations. 

Both the weak and the strong coupling approximations 
are based on a functional integral or field-theoretic rep- 
resentation [20| of the grand canonical partition function 
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of a system composed of fixed surface charges with in- 
tervening mobile counterions, and depend on the value 
of a single dimensionless coupling parameter S [l^ . The 
distance at which two unit charges interact with ther- 
mal energy fc^T is known as the Bjerrum length Is — 
eQ/ATreeoksT (in water at room temperature, the value 
is Ib ~ 0-7 nm). If the charge of the counterions is q then 
the aforementioned distance scales as q^ls- Similarly, the 
distance at which a counterion interacts with the surface 
charge a with an energy equal to fc^T is called Gouy- 
Chapman length, defined as /j, = eo/2TTqlB<J. A compe- 
tition between ion-ion and ion-surface interaction can be 
quantitatively measured with a ratio of both character- 
istic lengths S = g^Zs/^ = 27rg'^Zgcr/eo, which is known 
as the (Netz-Moreira) coupling parameter (l8l |. 

The meaning of this coupling parameter can be easily 
understood by considering the mean distance between 
counterions in the layer next to the charged surface. 
Assuming that electroneutrality is achieved within one 
Gouy-Chapman layer /x, the volume available per ion is 
Aira^ /3 = qeon/a and the mean distance between coun- 
terions a = /i(3S/2)^/'^. It follows that in the weak cou- 
pling case, defined by S <C 1, the width of the layer /i is 
much larger than the separation between two neighbour- 
ing counterions and thus the counterion layer behaves 
basically as a 3D gas. Each counterion in this case in- 
teracts with many others and the collective mean-field 
approach of the Poisson-Boltzmann type is completely 
justified. On the other hand in the case of the strong 
coupling limit, defined by S 3> 1, the mean distance be- 
tween counterions, a, is much larger than the layer width, 
meaning that the counterion layer behaves as a 2D gas 
[ri |. In this case the mean-field approach breaks down, 
each counterion moving almost independently from the 
others along the direction perpendicular to the wall and 
the collective effects that enable a mean-field description 
are thus absent. The two limits are characterised by a 
low/high valency of the counterions and/or a small/large 
value of the surface charge density. The range of validity 
of both limits has been explored thoroughly in [l^ . 
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FIG. 1: Schematic view on the charged surface with counteri- 
ons in the case of weak (lefthand drawing - H ^ 1) and strong 
(righthand drawing - H 3> 1) coupling hmit. \i is the Gouy 
- Chapman length and ax is the average separation between 
counterions. 

Formally the weak coupling limit can be straightfor- 
wardly identified with the saddle-point approximation of 
the field theoretic representation of the grand canonical 
partition function [20| , and is reduced to the PB theory in 
the lowest order. The strong coupling approximation has 
no PB-like correlates [l^ since it is formally equivalent 



to a single particle description, and corresponds to two 
lowest order terms in the virial expansion of the grand 
canonical partition function. The consequences and the 
formalism of these two limits of the Coulomb fluid de- 
scription have been explored widely and in detail (for 
reviews, see [HI, Q)- 

In this paper we thus embark on a study of interactions 
between macromolecules with pronounced helical motifs 
of the fixed charges that would allow us to formulate 
their interactions on the strict level of weak and strong 
coupling approximations. Our goal here is not to go into 
detailed modeling of the charge distribution as was done 
in the case of DNA by Kornyshev and Leikin [1] , neither 
do we intend to study the details of ss-DNA wrapping 
geometry and energetics in the case of decorated CNTs 
as was pursued by Lustig et al. Q. Our goal is more 
modest: we intend to asses the consequences of the weak 
and strong coupling dichotomy as it transpires through 
the interactions between macromolecules with a helical 
charge motif in the presence of counterions, and conse- 
quently to explore the specific features the helical charge 
motif adds to the interactions in both limits, studied ex- 
tensively heretofore TlJ. A similar, yet in fundamental 
respect different point of view, was recently advanced by 
Lee [H. 



II. THE MODEL 

We consider a model of two identical infinite parallel 
charged cylinders with radius a and interaxial separation 
R. The charge on both cylinders is distributed on a single 
helix with a pitch H. In order to avoid spurious diver- 
gences when solving the PB equation, we represent this 
helical charge distribution as a helical stripe with small 
but finite thickness. This can be justified by the finite 
size of the fixed charge groups on the macromolecular 
surface, as well as the finite size of the counterions. 

Both cylinders can be rotated around their axes by an- 
gles ipi and ip2. Without loss of generality we will put 
ipi = and ip2 = ^Po, which just shifts the origin of the 
coordinate-system. For continuous helices the rotation 
of the second cylinder by an angle tpQ corresponds to a 
translation by Az = {H/2Tr)(pQ along the direction of 
the axis of the cylinder. The charge of both cylinders is 
compensated by freely mobile counterions with valency 
q bathing the cylinders. In counterion-only case we ne- 
glect all colons. This approximation is relevant for low 
salt concentrations where the Debye screening length is 
much larger than the scales of interest. We consider the 
cylinders as impenetrable to counterions. In our model 
the cylinders are considered as hollow, so that they con- 
tain water and therefore have the same dielectric constant 
as the bulk water, e — 80. This can be justified for hol- 
low CNTs, whereas it is only a rough approximation for 
a DNA molecule, because it neglects image contributions 
Q. However, even if we include exact image effects this 
would still lead to a very rough approximation of a real 
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FIG. 2: Two cylinders with (single) helical charge distribution 
in a bath of counterions (not shown). 



DNA molecule because of unknown surface and satura- 
tion effects on the value of the local dielectric constant. 
We will not delve into these complicated and poorly un- 
derstood effects in this work. 

The amount of charge on both cylinders can be ex- 
pressed with a dimensionless Manning parameter (2ll.[23|. 
defined as 



Q = 



(1) 



where d is the longitudinal spacing between equivalent 
elementary charges along the cylinder and / b is Bjerrum 
length. Expressing the Manning parameter with mean 
surface charge density a, one remains with 
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(2) 



Note that at fixed surface charge the Manning parameter 
Q is proportional to counterion valency q. 

We will solve the statistical mechanics of the above 
model in two well defined limits [3]: the weak-coupling 
(WC) a.k.a. the Poisson-Boltzmann (PB) or the mean- 
field limit and the strong-coupling (SC) limit. Their per- 
tinent range of validity have been analysed thoroughly, 
see e.g. [14l | . These two coupling regimes delimit the ex- 
act thermodynamic properties of the system and thus de- 
scribe the extreme limits of its behavior in the parameter 
space. We will not include in our analysis the contribu- 
tion of the second order fluctuations around the mean- 
field solution [2^, as can be done for planar systems, 
since it is a lot more complicated in cylindrical geometry 
[23 |. while remaining numerically overall small. 



III. WEAK COUPLING LIMIT (MEAN-FIELD) 

The mean-field approach, also known as weak-coupling 
limit (S <C 1), is based on the Poisson-Boltzmann equa- 
tion for charged mobile counterions in solution [3] and is 



valid for S ^ I. This theory predicts that counterions 
will distribute in the space surrounding the macromolec- 
ular charge in accordance with the Boltzmann statistics 
that leads to the following equation for the mean electro- 
static potential 



-f3eoq4> 



eeo 



(3) 



where the constant A can be determined by the charge 
neutrality condition. Inside the cylinders the right-hand 
side of the above equation is 0. Using the dimensionless 
electrostatic potential, u = PeQq4>, we obtain a set of 
equations valid inside and outside the cylinders 



V^u = -Ce- 
V^u = 



(outside), 
(inside). 



(4) 



The boundary conditions on the surface in the presence 
of the surface charge density a* , that can exhibit spatial 
variation along the surface, are formulated in the stan- 
dard way via the normal component of the electric field 
strengths as 

where e is the (static) dielectric constants outside as well 
as inside the cylinders. In dimensionless variables this 
amounts to 



dr 



du 
dr 
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charged stripe, 
otherwise. 



(5) 



In order to avoid numerical divergences in the deriva- 
tives if the surface charge density a* is represented by an 
infinitely thin line, we represent the helical charge dis- 
tribution as a helical stripe with a finite thickness. As 
already noted this can be justified by realizing that a 
typical ion radius would be around 0.2 nm and therefore 
we assume the stripe of thickness 0.4 nm. 

Numerically we solve the PB equation in a finite 
bounding box geometry with dimensions L^xLyX Lz (see 
[2]) . The height of the box is always taken to be equal to 
one pitch, = H and the periodic boundary conditions 
are applied along the z direction to mimic infinitely long 
cylinders. Results depend on lateral box sizes Lx and 
Ly and should converge when these two sizes are large 
enough. 

Solving for the dimensionless potential u we can calcu- 
late the force acting between both cylinders via the stress 
tensor p, composed of the Maxwell part and the ideal gas 
(van't Hoff) part as 



Pik = ££0 EiEk 
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(6) 



The van't Hoff pressure of the counterions, p, equals their 
ideal gas pressure. The force between cylinders can be 
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evaluated by integrating the stress tensor over any closed 
surface embedding one of the cylinders: 



PikdS, 



k ■ 



(7) 



We choose to integrate over the mid-plane between the 
cylinders and over the outer edges of the box, as was done 
in With large enough simulation box the integral 

over the outer edge should become negligible. 

As stressed already in the introduction the helical 
charge model can be considered as an approximate ren- 
dition of either the ds-DNA or the ss-DNA/CNT hy- 
brid system. For the ds-DNA the structural parameters 
are: cylinder diameter of a = 1 nm and helical pitch of 
H = 3.4 nm, so that the Manning parameter is Q = 4.1q 
where q is the counterion valency. In the case of ss- 
DNA/CNT hybrid complex, the cylinder radius can typ- 
ically be in the range of 1.5 — 3 nm and the pitch from 
about 3.4 nm on 0,0]. Depending on these numerical 
values, the Manning parameter in the ss-DNA/CNT case 
is in the range of 2 — 10 for monovalent counterions. If 
not stated otherwise the counterion radius is taken to be 
i?c = 0.2 nm. Since ds-DNA case is in the middle of both 
ss-DNA/CNT extremes we will mostly delimit ourselves 
to DNA parameters as typical values used in the numer- 
ical calculations. They are used for illustrative purposes 
and not for an explicit description of ds-DNA per se. 

The coupling parameter for the ds-DNA case is S = 
2.8 q'^ in the presence of counterions with valency q. The 
results for WC as well as SC approximations do not de- 
pend directly on S but this parameter determines the 
range of validity of these approximate theories [l^ . For 
monovalent counterions the coupling parameter has a 
value of S = 2.8 which is small enough to justify the 
WC approach ^2^. But this approach certainly breaks 
down for tri- and tetra-valent case where S = 76 and 
180, respectively. 



A. Results 

In order to solve numerically the system of partial dif- 
ferential equations with spatially varying boundary con- 
ditions, Eqs. H and O we used the COMSOL Muhi- 
physics package, a finite element analysis and solver soft- 
ware, together with the Cholesky decomposition method 
for evaluating results on finite element mesh. 

Before we delve into helical charges we take a quick 
glance at the behaviour of two homogeneously charged 
cylinders [21]. First of all we want to explore the ef- 
fect of the finite bounding box in our numerical calcu- 
lations. The bounding box has a square base, Lx = 



L and height Lz 



H with periodic bound- 



ary conditions applied in the vertical direction. It is 
well known, that the counterion density around a sin- 
gle uniformly charged cylinder behaves asymptotically as 
r-2[l -f (Q - l)lnr/a]-2 [13, [11]. Therefore, we expect 
to see the same asymptotic behaviour for helical charge 



distribution at large distances also in our numerical cal- 
culations, with the proviso that Q — *■ 2Q for two cylin- 
ders. 

The density profile in our case decays very slowly, 
causing weak convergence with increasing box size. For 
this reason the finite-bounding box effects are important. 
This size dependence can be discerned from Fig [31 By 
constraining the space around the macroions we artifi- 
cially heap surrounding counterions into the box that 
should otherwise extend further away. At small inter- 




FIG. 3: Force per unit length - separation dependence for 
uniformly charged cylinders with ds-DNA parameters [q — 1) 



for different box sizes Lx 



L. The finite box size 



effect is of course more pronounced at larger separations. The 
inset directly shows the variation of the force with increasing 
box size at the distance R — 4 nm. One should note the 
logarithmically slow convergence of the result. 

axial separations R the force is obviously not as sensitive 
to the bounding box size as at larger separations. As 
shown in the inset of Fig [3l the convergence of the result 
is very weak when exceeds w 20 a. In all of the cal- 
culations presented below we used L = 25 a and we have 
to be aware that this result is slightly but not fundamen- 
tally different from the strict formal limit of an infinite 
box size. 

On the mean-field level the force between charged he- 
lices is always repulsive and decays with increasing inter- 
axial distance. In a symmetric system on the mean-field 
level it is well known that interactions are always repul- 
sive, a result that can be formulated as a general theorem 

M- 

On Fig 13] and [H] we present the results of mmrerical cal- 
culations of the interaction force between two cylinders 
with a helical charge motif as a function of the interax- 
ial separation as well as the azimuthal angle of mutual 
orientation. Increasing the Manning parameter (corre- 
sponding to higher counterion valencies) decreases repul- 
sive interaction but can not lead to a change of sign of 
the interaction. The maximal angular variation at a sep- 
aration of i? = 2.3 nm amounts to only 16% of the total 
force per unit length, and is a lot smaller if the cylinders 
are even further apart. 
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R [nm] 



FIG. 4: Force per unit length - separation dependence for two 
lielices (ds-DNA parameters) at various mutual azimuthal an- 
gles ifio and valencies. The corresponding Manning parameter 
for 5 = 1 is Q = 4.1 and for q = 2 it is Q = 8.2. 



The azimuthal angle dependence of the force for two 
helical charge distributions, Fig O thus shows only a 
slight variation, being larger at smaller interaxial sep- 
arations. This is a completely general conclusion of the 
WC analysis. The azimuthal modulation of the force has 
a maximum at a configuration corresponding to fo = tt 
and a minimum at (/Sq = 0. It also shows mirror sym- 
metry along 1^0=71 line, so that configurations ipo and 
2TT — ipo have exactly the same interaction. 
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FIG. 5: Azimuthal angle dependence of the force per unit 
length at different interaxial distances for ds-DNA parameters 
and monovalent counterions (Q = 4.1). 



IV. STRONG COUPLING LIMIT 

In the limit of S ^ 1 the mean-field ansatz breaks 
down and a different kind of approach is needed. Here, 
virial expansion of the partition function in terms of fu- 



gacity A' can be used [l^l , leading to the strong-coupling 
limit. For our purposes only the expansion of partition 
function Zo to the first order will be used, that effectively 
amounts to 



7 _ 7(0) 

— 



0(A'2) 



(8) 



Here the zeroth-order term corresponds to the bare elec- 
trostatic interaction energy of charged cylinders, and the 
first-order term corresponds to the one particle counte- 
rion contribution to the partition function. For details 
of the strong-coupling approach, see (T3 | and references 
therein. 

To obtain the free e nerg y in the SC limit we follow the 
same procedure as in 29]. The first step is to evaluate 
the grand potential from the grand partition function and 
then by using a Legendre transformation we get the free 
energy per counterion as 



f3T 



N 



In 



(9) 



The first term in the above expression, Wq, represents 
the bare electrostatic energy between charged helices. It 
can be evaluated simply by the integration of charge con- 
tributions on both helices with the unscreened Coulomb 
kernel 



N 



2H 



" dzdz' 
s{z,z')' 



(10) 



Here, s{z, z') represents the distance between two points 
on different helices with coordinates z and z', respec- 
tively. Helices are represented as (infinitely) thin lines 
and not as a stripes, as was the case in the WC case. 
On the SC level the counterion hard core radius will be 
taken into account via its excluded volume in the parti- 
tion function [2^ . The second term in Eq. [9] is a one- 
particle contribution to the free energy. Here /3u is the 
electrostatic energy of a single counterion with charge 
eoQ in the presence of two charged helices. It is obtained 
simply by integration of the Coulomb kernel along both 
helices as 



1 



1 



Sl(z') S2(Z') 



dz\ 



(11) 



where si and S2 are distances from the counterion- 
coordinate to the point on helix 1 and helix 2 defined 
as 

si,2(p,</?, z;^') = (12) 

\J a? + — 2ap cos(fcz' + (^1,2 — </?) + (z — z'Y 

The angles (pi and Lp2 represent orientations of both cylin- 
ders around their symmetry axes. Here, again as in PB 
case, the convention Lpi = Q and = (po is used without 
loss of generality. 
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The integral in Eq. [9] should be in principle eval- 
uated over the entire space available to the counteri- 
ons, but we again delimit integration inside the confin- 
ing box. The height of the box is again taken as one 
helical pitch, — H, due to the periodicity of the po- 
tential u{r) in z-dircction. Any multiple value of H in 
Lz would produce only an additive constant in the final 
free energy result. The problem of finite integration box 
in the strong-coupling limit is not so critical as in the 
Poisson-Boltzmann case, since counterion density decays 
very rapidly with radial distance. In the case of a single 
uniformly charged cylinder the density decays as r~^^ , 
so the asymptotic behaviour at large distances for two 
cylinders with a helical charge motif behaves as r^*'^ . 
Typically, Q > 4 for DNA so the volume integral con- 
verges very rapidly. A rectangular box of sides = 7a 
and Ly — 3a was used, which is large enough to produce 
results very close to the infinite box size limit. 

Hard core repulsion of cylinders as well as finite 
counterion-size are taken into account via the hard-core 
radius of the cylinders, equal to a -I- Re- Here Rc is coun- 
terion radius and a the bare radius of the cylinder. 

The interaction force between both cylinders is ob- 
tained by taking the derivative of the free energy with 
respect to the interaxial separation 



F 



'dR' 



(13) 



Obviously, the force scales with the length of the cylin- 
ders, therefore we express the force per length of cylin- 
ders, F/L. 
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FIG. 6: SC force per unit length - separation curves for vari- 
ous orientations tpo in the case of ds-DNA parameters of a = 1 
nm and H = 3.4 nm in the presence of monovalent counteri- 
ons. The corresponding Manning parameter is Q — 4.1 and 
the coupling parameter H ~ 3. The inset shows force per unit 
length - separation curves for orientation (po — tt for different 
cylinder radii, corresponding to the numerical values of the 
parameters of the ss-DNA/CNT hybrids. Qualitatively the 
interactions are the same. 



Similar to the case of interacting planar surfaces, for 
very small interaxial spacings the counterions accumu- 
late in the space between the two cylinders, creating a 
repulsive component to the total interaction due to their 
osmotic pressure. The osmotic contribution can be larger 
than the electrostatic correlation attraction and a net re- 
pulsion ensues for the total interaction between helices 
at very small separations. 



A. Results 

The interaction force in the case of the strong-coupling 
limit shows a much richer behavior than in the weak- 
coupling (PB) case. As seen from Fig [S] the interaction 
force for various parameters, including the ds-DNA as 
well as ss-DNA/CNT values, exhibits similar qualitative 
behaviour. Again we will use ds-DNA parameters for 
illustrating our analysis. We present the SC results for 
monovalent as well as polyvalent counterions even though 
small valencies q that correspond to small coupling pa- 
rameter S are not relevant for the SC theory. The aim 
here is to illustrate different behaviour of the interac- 
tion for counterions of different valencies, i.e. of different 
Manning parameters. 

We find three different regimes in the behavior of the 
interaction as a function of interaxial separation, Fig [6l 
At large separations the SC force is always attractive. 
This is due to the fact that each counterion is localized 
at one or the other cylinder, effectively creating a cor- 
relation hole, causing net correlation attraction between 
the helices. This is indeed very similar to the interaction 
of planar surfaces [l3|- At larger interaxial spacings the 
SC theory breaks down and the WC repulsion takes over. 
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FIG. 7: The hard core counterion-radius effect. Depletion 
region, where the force is reduced, appears at distance smaller 
than counterion diameter. 

The highest value of osmotic pressure is reached when 
the separation between cylinders becomes equal to the 
counterion diameter 2Rc, see Fig [T] This leads conse- 
quently to the highest repulsive force between the cylin- 
ders. At yet smaller separations, the depletion effect sets 
in, in the sense that the counterions can not penetrate 
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FIG. 9: Force-separation curves for various counterion valencies q at different orientations ipo (ds-DNA parameters). At large 
distances the force is independent of the valency. 



anymore the inter- cylinder space, and are thus depleted 
from the spatial region between the two apposed cylin- 
ders, diminishing their osmotic pressure in that region. 
Thus, the counterion contribution to the osmotic pressure 
is reduced, i.e. the osmotic repulsion as well as counte- 
rion electrostatic interactions are both reduced. In all 
these situations the bare electrostatic repulsion between 
the two helices remains unaffected. Because of this, at 
a separation R* where the net force is zero F — 0, the 
bound state between the cylinders occurs. 

The relative orientation of the two cylinders plays a 
crucial role in determining the nature of the inter-helical 
interaction. The force at orientation ipo = n now appears 
an order of magnitude larger than at orientations or 
7r/2. This is quite different than in the WC case, where 
the analogous variation was much smaller, and indeed 
almost negligible. The case of two uniformly charged 
cylinders is found somewhere in between both extremal 
behaviors, see Fig[Sl 
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FIG. 8: Force per unit length - angle curves for various inter- 
axial distances. DNA parameters for q = 1 are taken. 

From Fig [8] we also deduce that the extremum of the 
force is reached at (po = tt. This could be explained 



by the fact that at ipQ — tt the charges on helices are 
apposed and can approach at a smallest local separa- 
tion. The electrostatic potential between the cylinders is 
therefore higher than for other values of ipo , which causes 
stronger localization of counterions in between. Stronger 
localization in its turn leads to larger attraction at large 
distances as well as greater osmotic repulsion at smaller 
distances. 

As seen from Fig [9] the counterion valency, i.e. Man- 
ning parameter Q, can also have quite a dramatic effect 
in the SC limit. Note that increasing the valency q also 
increases the Manning parameter in the same propor- 
tion. At large distances the force is independent of the 
valency. This makes sense, since at large distances each 
counterion is located either at one or the other cylinder. 
In this case counterions cannot sterically interact with 
both cylinders, therefore the osmotic repulsive contribu- 
tion vanishes. Only electrostatic interaction thus remains 
that depends on the net counterion charge but not on the 
amount of counterions, i.e. it does not depend on coun- 
terion valency. Note, that higher valency means smaller 
amount of counterions since the net counterion charge 
remains the same due to electroneutrality. 

The influence of counterion valency is significant also 
at smaller interaxial separations where osmotic compo- 
nent plays an important role. Each counterion con- 
tributes to the osmotic pressure in the vicinity of cylin- 
ders. Its contribution is larger for higher electrostatic 
potential around cylinders as well as larger valency. 
Namely, counterions with larger valencies are more at- 
tracted and localized in the vicinity of cylinders. This 
dependence of the osmotic force on the valency is non- 
trivially connected with the electrostatic potential. The 
net osmotic force in SC is the product of singlc-counterion 
contribution and the amount of counterions N. Due to 
electroneutrality condition, this amount is then inversely 
proportional to the valency, i.e. N oc q^^ . The net os- 
motic contribution can thus be an increasing as well as a 
decreasing function of q. 
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Fig[U]shows the dependence of the force per unit length 
on the interaxial spacing for various orientations as well 
as various counterion valencies. In (a) the repulsive os- 
motic contribution ai ipo ~ is apparently higher for 
smaller valencies and its value per counterion increases 
slower than q. We found out, that such a trend also 
appears if the cylinders are uniformly charged, the case 
analyzed in 21]. In the (po = n case, the opposite is 
true (Fig [5] (c)). In this case, the segments of charged 
helices approach to the smallest mutual separation, caus- 
ing much higher electrostatic potential in between, that 
enhances counterion localization which furthermore leads 
to higher osmotic repulsion for larger valencies. Thus on 
increasing q the enhanced localization leads to more pro- 
nounced growth of the osmotic contribution per counte- 
rion so that the net osmotic force increases with increas- 
ing q. In both configurations ipo = and ipo ~ ir the 
largest potential appears exactly in the middle of both 
cylinders which induces a maximal counterion density at 
the same position. This is however not the case for other 
orientations. 

At the configuration (pQ — tt/2 (Fig [5] (b)) an even 
more dramatic phenomenon occurs. Here, the highest 
potential appears slightly out of the interaxial plane. The 
counterions that are localized slightly away from the cen- 
tral region between the helices also contribute less to the 
osmotic pressure. For q = 1 {Q = 4.1) case when the 
amount of counterions is still large enough and localiza- 
tion not so strong the osmotic pressure prevails over the 
electrostatic attraction at smaller separations. But for 
g > 2 (Q > 8) the attractive electrostatic interactions 
prevail at all separations, leading to a collapse transition 
of both cylinders into a bound state, corresponding to the 
closest approach distance, which is the double radius of 
the cylinders 2a. As seen from FiglTUb. the bound-state 
distance R* changes only slightly for small Manning pa- 
rameters, corresponding to monovalent counterions for 
the ds-DNA case. The occurrence of this transition can 
be understood from the free energy-interaxial distance 
plot in Fig [To] (b) . The free energy exhibits two min- 
ima at small orientations <po, with the minimum at finite 
separation being the global one. On increasing ipo the 
free energy minimum located at the closest separation 
R = 2a becomes deeper, leading to a discontinuous tran- 
sition into a bound state. At even larger orientations 
the minimum continuously moves away from the closest 
separation to finite separation. 

Around iy9o = and ipo — n the bound-state distance 
is very close to the case of two uniformly charged cylin- 
ders. There, a simple analytical expression can be de- 
rived, R* — 2a ^ 2Rc + |/i [30|. The difference in separa- 
tions between helical and uniform charge distribution is 
within 20% for monovalent and less than 10% for poly- 
valent case for aforementioned orientations. But it com- 
pletely fails for intermediate (po for polyvalent case. 

We also want to briefly comment on the behaviour 
for small Manning parameters. According to Naji et al. 
[30| . a significant dependence upon the size of the con- 
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FIG. 10: Collapse transition into a bound state for ds-DNA 
parameters, (a) Bound state interaxial separation as a func- 
tion of mutual orientation tpo. For polyvalent ions the rapid 
jump appears at certain angle, corresponding to phase transi- 
tion, (b) Corresponding free energy-distance plot for divalent 
counterions {q = 2) at different orientations. Two minima (for 
the first three cases) are clearly visible which are responsible 
for the phase transition. 



fining box appears for Manning parameters smaller than 
Qc = 2/3 in the case of two uniformly charged cylin- 
ders. This occurs as a result of the dilution of the coun- 
terion cloud around weakly-charged cylinders. Below the 
critical value of the Manning parameter, Qc, the equi- 
librium distance scales with the box size approximately 
as i? ^ Lx/y/TT. We found that helical charge distribu- 
tion in our case exhibits the same behaviour below the 
critical Manning parameter Qc = 2/3. This is expected, 
since the charge distribution details are not important 
at large separations. Such small values of Q do not of 
course correspond to DNA parameters. 

Within our model there is no attraction between the 
helices on the mean-field level, so there can not be any 
bound state. In the ds-DNA model of Kornyshev and 
Leikin [1, [3l| the counterions are assumed to be local- 
ized and fixed on the surface of the helix, leading to a 
charge separation between negative charges of the phos- 



9 



phates and positive charges of the (fixed) counterions. 
In that case of course, even on the Hnearized PB level, 
one can still see attractions due to charge separation, 
leading to a bound state with a typical equilibrium sep- 
aration of 0.5 nm, depending on the parameters - close 
to the experimentally reported value in bundles of DNA 
molecules [slj. Fixing the counterions on the DNA sur- 
face of course assumes an additional specific binding in- 
teraction, not present in our model where the only in- 
teraction is Coulombic, so the results for the Kornyshev- 
Leikin and our model can not be compared directly at 
this point. 

If the helices are rotationally unconstrained, they will 
always settle into an orientation that corresponds to a 
minimum of their free energy. In the weak-coupling case 
this orientation always corresponds to (fo = 0, since in 
this configuration the helical segments are locally maxi- 
mally apart and the (repulsive) free energy is thus min- 
imal. This was already emphasized by Kornyshev and 
Leikin [1, [3^ in the linearized PB approach for two sin- 
gle helices. 

However, the behavior of this optimal azimuthal angle 
in the SC case is very different and in many respects op- 
posite to that observed on the WC linearized PB level 
[1, [s^]. At large distances the optimal angle is tpo — tt. 
This is in fact the most unfavorable configuration in the 
WC case. Correlation effects in SC dramatically change 
the behaviour of the system. Here, counterions mediate 
electrostatic correlation attraction between both helices 
and thus they tend to orient themselves so that the lo- 
cal helical segments approach as much as possible. This 
happens exactly at the ipo = tt configuration. We must 
emphasize here that only the electrostatic (correlation) 
interaction plays a role here while the repulsive osmotic 
contribution to angular interaction is negligible since it 
cannot cause any torque, acting always perpendicular to 
the cylindrical surface. That is the reason that the op- 
timal angle remains ipo = tt even though the net force 
turns repulsive! 
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FIG. 11: The interaxial distance dependence of the optimal 
angle between two helices with DNA parameters. The critical 
separation where the rapid change of optimal angle occurs 
equals the counterion diameter, R — 2a = 2Rc = 0.4 nm. 



Everything said heretofore remains valid for larger 
surface-to-surface cylindrical separations, when com- 
pared to the counterion diameter 2i?c- At smaller sepa- 
rations the counterions are depleted away from the mid- 
region between the cylinders. The localization of counte- 
rions is then out-of-center and causes a torque on helices 
resulting in a different value of the optimal angle ipQ . As 
shown on Fig [TT] the optimal angle above the critical sep- 
aration {2Rc~2a = 0.4 nm) is tt. Below this separation a 
substantial variation of the angle can be discerned. Note, 
that configurations Lp^ and 27r — Lp^ are equivalent, so we 
always plot the one with the smaller azimuthal angle. 

The reason for the variation of the optimal angle in 
SC at smaller separations is the depletion effect that is 
not present in the mean-field case. Similar behaviour of 
the optimal angle below a critical value of the interaxial 
separation was also found on the linearized PB level for 
multistranded helices. There, the spontaneous symmetry 
break occurs at a critical separation where the optimal 
angle switches from (po = to higher values 



B. Regime of applicability of SC results 

So far, we evaluated numerical results for the strong 
coupling limit for charged helices without any justifi- 
cation of its validity. In a planar system with a finite 
value of the coupling parameter S, the asymptotic strong- 
coupling results hold exactly as long as the surface sep- 
aration 5 is smaller than the typical lateral distance be- 
tween counterions a± . This condition in fact yields a 
simple and generic criterion identifying the regime where 
strong-coupling attraction is expected to emerge between 
two like-charge planar macroions. It was originally sug- 
gested by Rouzina and Bloomfield [Ts'l. In the case of 
uniformly charged cylinder with a large coupling param- 
eter S, counterions tend to line up on opposing surfaces 
of the cylinders and along axes forming a correlated inter- 
locking pattern [34]. Typical distance between counteri- 
ons along the axis of the cylinder, may be estimated 
from the electroneutrality condition 



2A ■ 



(14) 



The strong coupling limit is expected to become valid 
when the surface-to-surface distance of the cylinders, 
6 — R — 2a, becomes smaller than the distance between 
counterions, i.e. 6 < az, leading to the criterion [s^ 



(15) 



where Qi is the Manning parameter for monovalent coun- 
terions. Note that this critical separation grows linearly 
with valency q. Taking the parameters for ds-DNA, with 
Qi = 4.1 and Bjerrum length Ib = 0.7 nm the values 
of critical distance 5 are in the range from 0.2 nm to 
0.7 nm when the valency q goes from from 1 to 4. This 
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calculation suggests that tri- and tetra-valent counteri- 
ons do justify the SC approach at separations around 0.5 
nm and below, where most of the interesting phenomena 
described above occur anyhow. But according to this cri- 
terion, monovalent and divalent counterions can not be 
described within the SC limit for the ds-DNA parameter 
case. 
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FIG. 12: Schematic representation of SC and WC regions 
for ds-DNA parameters. (ySo = in the presence of tetrava- 
lent counterions, g = 4. In the intermediate separations no 
approach is valid and one would have to take recourse to ex- 
tensive simulations. 

Note that we used the criterion for uniformly charged 
cylinders and that helical charge pattern might change 
this criterion since counterions are localized and line up 
differently. The regime of SC-validity might also depend 
on orientation of molecules. But so far we must rely on 
this simple and rough criterion. 

V. CONCLUSION 

In this paper we have analysed the electrostatic in- 
teractions between two cylinders with a helical charge 
distribution in the presence of counterions. We used the 
approach a la Netz and coworkers by explicitly consider- 
ing the limits of weak and strong coupling, defined by the 
value of the electrostatic coupling parameter S, which is 
assumed to be small in the former and large in the latter 
case. 

The weak coupling case, or equivalently the mean-field 
case, is relevant when counterion valency is small, i.e. 
g = 1, or at large inter- helical separation. The distribu- 
tion of counterions and the corresponding electrostatic 
potential is governed by the Poisson-Boltzmann equa- 
tion. The force between two charged helices has been 
evaluated by the stress tensor method and appears to be 
overall repulsive in the WC case, a result known from pre- 
vious works [ll, [3, HI] . Helical charge distribution here 
makes up for only a slight angle modulation with compar- 
ison to the uniformly charged case. The largest repulsive 
force appears at (fo = tt mutual cylindrical orientation 



and the smallest at Lpo = 0. Increasing the counterion 
valencies q tends to decrease the amount of counterions 
in the system and the corresponding osmotic contribu- 
tion as well as the net interaction is therefore smaller 
but it never changes the sign. 

In contrast to weak coupling, the strong coupling limit 
exhibits a much richer and a lot more interesting behav- 
ior. The SC theory is effectively a one-particle theory 
[3] that neglects all interactions among counterions and 
only takes the interactions between the counterions and 
the fixed charges on the macromolecular surface, e.g. the 
ds-DNA surface, into account. This approach is justified 
if the distances between counterions are much larger than 
their distances between the cylinders. In other words, 
SC limit is valid for large counterion valencies q and for 
small separations between the cylinders. Particularly for 
ds-DNA parameters, the SC approach for q = 3 and 4 
would be valid for surface-to-surface separations under 
~ 0.5 nm which is comparable to typical counterion di- 
ameters. 

As we have been arguing above, the SC force be- 
tween two cylinders is always attractive at larger sepa- 
rations, but can become repulsive at smaller separations 
due to the osmotic counterion contribution. At sepa- 
rations smaller than counterion diameter the counterions 
are depleted away from inter-helical region which reduces 
the osmotic repulsion. At the separation where the net 
force is zero, a bound state occurs. This bound-state sep- 
aration approximately equals to the counterion diameter, 
typically 0.4 nm. 

On the SC level the relative orientation of cylinders has 
a drastic effect on interactions between cylinders with 
a helical charge distribution. In many respects the in- 
teractions between helices in this limit are opposite to 
those observed on the WC linearized PB limit [o, '33^ . At 
the relative orientation tpQ = tt the helices on different 
cylinders approach locally to smallest possible separation 
causing high electrostatic fields that attract counterions 
more than in other configurations. The interaction in 
the orientation ipQ = tt is therefore an order of magni- 
tude larger than for tpo = ot ipQ = tt / 2 at small separa- 
tions. At large separations the interaction is independent 
of orientation and specific charge distribution. There it 
expectantly approaches the uniformly-chargcd-cylinders 
results. 

At well defined conditions the osmotic counterion con- 
tribution can be reduced to such an extent that the 
electrostatic correlation attraction prevails and pulls the 
cylinders into a bound state at their closest possible sep- 
aration, R = 2a. Particular for ds-DNA parameters this 
happens for multivalent counterions, q > 2, and in de- 
fined orientation range (fo O.Stt — O.Gtt. The counte- 
rions at these conditions are localized away from the in- 
teraxial plane and the osmotic contribution is therefore 
sufficiently reduced for this to happen. 

If two helices can freely rotate around their long axes 
they will settle in the orientation that minimizes their 
free energy. At large distances the optimal angle for SC 
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is ipo = TT, which is the most unfavorable configuration 
in the WC case. Counterions that mediate electrostatic 
attraction between both helices tend to orient helices so 
that the local helical segments approach as much as pos- 
sible. When surface-to-surface separation is smaller than 
the countcrion diameter, the depletion effect localizes 
counterions out-of-center and the optimal angle changes. 

Our approach in many respects complements the anal- 
ysis performed by Kornyshev and Leikin Q, which is 
based on the assumption of complete and ion-specific ad- 
sorption of polyvalent counterions onto the surfaces of the 
helical molecules. Rather, in the case investigated here, 
the only interaction between the counterions and the 
macromolecular surface is electrostatic in nature. The 
fact that even in this case there exists a strong angular 
dependent attraction between two cylindrical molecules 
bearing a helical charge motif is quite remarkable in itself, 
and results of our analysis will certainly be relevant in the 



case of helical charge distributions where there is no rea- 
son to expect any specific non-electrostatic interactions 
between the counterions and the macromolecular surface. 
Further work is in progress with the aim of understand- 
ing the effects of finite concentration of monovalent salt, 
apart from the presence of polyvalent counterions, on the 
behavior of this system. 
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